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Abstract 

Various quantities of an attractively interacting fermion system at the unitary limit are deter- 
mined by extrapolating Monte Carlo results of low-density neutron matter. Smooth extrapolation 
in terms of l/(/ci?ao) {kp is the Fermi momentum, and oq is the "^Sq scattering length) is found 
with the quantities examined: the ground-state energy, the pairing gap at T ~ 0, and the critical 
temperature of the normal-to-superfluid phase transition. We emphasize proximity of the physics 
of low-density neutron matter to that at the unitary limit. The extrapolated quantities are in a 
reasonable agreement with those in the literature. 
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I. INTRODUCTION 

In our previous paper [])], we reported a Monte Carlo calculation of thermodynamic 
properties of low-density neutron matter by using nuclear effective field theory (EFT) j^. 
As pairing in neutron matter is strong, neutron matter is a strongly correlated fermionic 
system. We firmly established that low-density neutron matter is in the state of BCS- 
Bose-Einstein condensation (EEC) crossover [3] instead of a BCS-like state, the standard 



description in nuclear physics {4]. 

The crossover state of low-density neutron matter is actually an expected one from studies 
of the BCS-BEC crossover over the past decade. The pairing strength of an attractively 
interacting (therefore unstable) fermion system is characterized by a product of the two 
physical parameters, the Fermi momentum kp and the (S-wave) scattering length [5]. In 
terms of = l/(A;^ao), the state of the BCS limit is realized at rj —00, that of the BEC 
limit at ?7 +00, and that of the unitary limit at 77 ^ 0. (We use the sign convention < 
for fermions attractively interacting with no bound state.) A fermion pair at the limit forms 
a zero-energy bound state, thereby yielding an infinitely long scattering length. Because the 
infinitely long scattering length generates no classical scale, the fermion system composed 
of the pairs is expected to have a universal feature. In recent years, much attention has 
been paid to the physics at the unitary limit in the fields of both condensed matter and 
atomic physics; for example, Refs. 0,0], and references therein. The ^5*0 neutron-neutron 
scattering length is negative and large in magnitude, oq ~ —16.45 fm Q]. For a density of 
about (10~^ — 10~^)po (Po = 0.16 fm~^ is the nuclear matter density), rj takes the value of 

-0.8<r/<-0.2. (1) 

Furthermore, the neutron matter of the above density is well described by the EFT lattice 
Hamiltonian that is identical to the attractive Hubbard model jo]. Equation ([1]) thus suggests 
that the physics of low-density neutron matter would be similar to that of the unitary limit. 
Note that the similarity is expected only in the limited range of the density. For 77 <^ —0.8 
(?7 — >■ —00), neutron matter would be in more of a BCS state in this quite low density. 
For rj > —0.2 [t] 0), the physics of neutron matter is much affected by the next-order 
(repulsive) term in the EFT potential, and deviates from that of the unitary limit. 

Physically, neutron matter never reaches the unitary limit 77 = 0. But our previous 
Monte Carlo results at the leading-order (LO) ^ are based on the Hubbard model with 



the parameter values chosen (translated from kp and oq) to be suitable for neutron matter. 
A Monte Carlo calculation of a fermion system at the unitary limit could be carried out 
by repeating our neutron-matter calculation but with the par ameter values adjusted to 



the limit. This approach has been used by Lee [10|, [111, [l^, [ISj, based on the Lagrangian 
formalism. While it may appear to be a straight-forward application, the computation is 
quite time consuming. Instead, by exploiting the proximity of low-density neutron matter 
to the unitary limit, we follow here another procedure originally proposed in Ref. [l^, an 
extrapolation to the unitary limit from low-density neutron matter by use of our previous 
Monte Carlo results at the LO. 

The underlying assumption in this procedure is that the physical quantities of interest 
are, at least numerically, slowly varying functions of rj and smoothly reachable to the unitary 
limit, as Eq. ([1]) suggests. As will be seen in Sees. Ill and IV, we find that this is indeed the 
case. The quantities examined are the ground-state energy Egg., the pairing gap at T 
A, and the critical temperature of the normal-to-superfiuid phase transition Tc, all of which 
are made to be dimensionless by taking ratios with ep = kp/{2M) (M is the neutron mass). 
Note that only kp provides a classical dimension at the unitary limit in the system. 

This paper is organized as follows: After the Introduction. Sec. II briefly summarizes our 
computational method, which will help keep the discussions in later sections coherent. In 
Sees. Ill, IV, and V, we describe the extrapolation procedure of Eg g , A, and T^, respectively. 
The description of Eg ^.jn Sec. Ill is somewhat detailed, since this quantity was not examined 



in our previous paper Our conclusion is found in Sec. VI. 



II. MONTE CARLO COMPUTATION FOR THE STANDARD PARAMETER 
SET 

The neutron-neutron [nn) interaction in the EFT Lagrangian includes all possible terms 
allowed by symmetries of the underlying theory of quantum chromodynamics (QCD) j^. 
The nn potential is in the momentum expansion form 

V{p', p) = co(A) + C2(A)(p2 + p") + ■■■- 2c2(A)p ■ p' + ■ • • , (2) 

where p and p' are the nn center-of-mass momenta, and A is the regularization scale. For a 
description of low-density neutron matter, we use a truncated potential including only the 
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first term in Eq. ([2]). The truncated EFT Hamiltonian on a three-dimensional cubic lattice 
H then takes the form of the three-dimensional attractive Hubbard-model Hamiltonian , 



= -t ^ cl^Cj^ + 6tY^ fii^ + ^co(a) ^ hi^fiii, 



(3) 



where a is the lattice spacing, t = l/(2Ma^) is the hopping parameter (M, the neutron 
mass), and denotes a restriction on the sum of all neighboring pairs. and Cjo- are the 
creation and annihilation operators of the neutron, respectively (cr =t, i), and hi^j = cj^Cia 
is the number operator with the spin a at the i site. 
The lattice spacing a is directly related to A as 



A ~ -. 

a 



(4) 



A should be set larger than the momentum scale, below which the truncated form of the 
potential is valid 15|]. A can be chosen to be smaller but should be at least 



A>p (5) 

for the momentum p at which the physics is studied [l^. In our case, we take p ~ kp, 
because the momentum scale corresponding to the excitation energy of interest is much less. 

The EFT parameter co(a) in the lattice regularization with the finite lattice spacing a is 
given as l3] 

47r/l 26^-' 
Co a =Jj( ' (6 

where 



9i = 



dx dy dz 



— (cos X + COS y + cos z) 



1.58796- 



(7) 



is a form of Watson's triple integral 18|, with p denoting the principal value of the integral. 



At the unitary limit (|ao| oo), we have 



co{a)/{aH) -An/9i = -7.91353 • • • , 



which agrees with the value found in the literature 6|. 



(8) 



In our Monte Carlo calculation at the LO 



we carried out an extensive Monte Carlo 



lattice calculation for the standard parameter set, which consists of the three values of the 



neutron matter density p with the lattice filhng (or the site-occupation fraction) n set to be 
1/4. The reasoning underlying this choice is somewhat involved, and we refer the reader to 
our previous paper jl|. Here, n is defined as 

n = o!^P = ^^{cl^Ci^}, (9) 

^ 1,(7 

where Ng is the lattice size (the total number of the three-dimensional lattice sites). Note 
that in this work, we specify the density of the interacting fermion system, p, using kp 
through the relation p = kp/lSTT"^). We emphasize that p here is an expectation value 
obtained by our Monte Carlo calculation, as Eq. iQ shows. 

The three densities of the standard parameter set are kp = 15, 30, and 60 MeV. The 
corresponding three values of t] are listed in Table I, together with those of a and co(a), which 
follow from the kp values with n = 1/4. Though not used in the Monte Carlo calculation, 
the value of co(a) at the unitary limit, Eq. ([8]), is also shown in the table for comparison. 
Note that a in the standard parameter set satisfies the EFT regularization condition, Eq. ([5]). 

TABLE I: Standard parameter set and Co{a)/{a^t) at the unitary limit rj = 0. 

rj Co /{ah) kp (MeV) a (fm) A (MeV) 
-0.7997 -5.308 15 25.64 24.18 
-0.3999 -6.354 30 12.82 48.36 
-0.1999 -7.049 60 6.409 96.73 

0.0000 -7.914 



The Monte Carlo calculations were performed on cubic lattices of Ng = 4^, 6'^, 8^, and 
10^ by the method of determinantal quantum Monte Carlo [2^, commonly used in 
condensed-matter physics. Using the data on the four different Ng^s, we apply the method 
of finite-size scaling to extrapolate to the thermodynamic limit. We take the continuum 
limit {n 0) by extrapolation using Monte Carlo data for various n's on the A^^ = 6^ 
lattice. In the following section we elaborate on how we determine Eg g^ at these limits and 
how we then extrapolate to rj —>■ 0. 
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III. GROUND-STATE ENERGY Er, 



Following common practice, we express the ground-state energy per particle of neutron 
matter Eg ^^ in terms of the energy parameter ^, 



E. 



g.s. 



5E, 



g.s. 



(10) 



EpG 3ei? 

where Epc is the ground-state energy per particle of the corresponding noninteracting sys- 
tem. ^ is expected to be of a universal character at the unitary limit. As done in our 
previous paper, we determine ^ at the three values of t] from the LO lattice calculations, 
first by taking the thermodynamic limit and second by applying the continuum limit. After 
carrying out the two steps, we extrapolate C, to r] = 0. 
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FIG. 1: Finite-size scaling of the energy parameter ^. The thermodynamic limit is at A^^ — > oo. 
The Monte Carlo data of with statistical uncertainties are shown at the Fermi momentum of 

— 1/2 

/ci? = 15 (triangle), 30 (square), and 60 MeV (circle). The dotted lines are the A'^s linear fits, 
Eq. (HH), to the data of iV^ = 4^ 6^ 8^, and 10^. 



First, we determine ^ of the thermodynamic limit at the three rj by applying the method 
of finite-size scaling using Monte Carlo data for the four lattice sizes, A^^ = 4^, 6^, 8^, and 
10^. The data used for the determination are shown with statistical uncertainties in Fig. [H 

As in the case of A discussed in Ref. jl|, the Ng dependence of ^ is found to be weak. The 
scaling exponent is difficult to determine, and the best fit to the Monte Carlo data results 



6 



in an exponent with a large uncertainty, essentially being indefinite. As shown in Fig. [T], 

— 1/2 

we find the choice of the Ng scaling power of ^, ~ iV," ^ (the same as that of A [1]) works 
reasonably well though not quite ideally. Because of the limited number of our Monte Carlo 

— 1/2 

data, we decided to proceed with the analysis using the A'^s scaling. With this scaling, 
the best fits are found to be 



-0.8, 



-0.4, iV. 



-0.2, 



-0.37(13) N-^/^ + 0.5784(35), 
-0.162(84) iV7^/2 ^ 0.4339(75), 
-0.131(63) ^"7^/2 ^ 0.3400(80), 



(11) 



and are shown with the Monte Carlo data in Fig. [H The last constant in each equation in 
Eq. ffTTj) is ^ in the thermodynamic limit (A^^ ~* oo). 

The preceding best fits are performed using the jackknife method (often used in the lattice 
QCD data analysis [2])]). The method is used to obtain all best fits in this work. 




FIG. 2: Energy parameter ^ as a function of the lattice filling n. Monte Carlo data (solid circles 
with statistical uncertainties) are for n = 1/16, 1/8, 3/16, 1/4, 3/8, and 1/2 on the A^^ = 6^ lattice 
for kp = 60 MeV. The dashed line, Eq. (jl2p . is the best fit to the data. 



Second, we take the continuum limit through n 
Ref. 



in a procedure similar to that in 



6|, as discussed in Ref. [l]. As shown in Fig. [2], the Monte Carlo data for the lattice 
fillings of n = 1/16, 1/8, 3/16, 1/4, 3/8, and 1/2 are best fit by 

C(n) = -0.132(24) n^/^ + 0.409(15). (12) 



Note that n = 1/2 corresponds to the quarter filhng of the lattice. The data are taken on a 
lattice of A^^ = 6^ at the Fermi momentum kp = QO MeV. Equation f|T2|) gives the ratio of 
^{n 0) and ^(n = 1/4) as 1.26(6). 

As discussed in Ref. jl|, such a ratio is expected and is also confirmed to depend weakly 
on Ng and kp in the case of A at T 0, T^, and the pairing temperature scale T* . Here, 
we assume that ^ also weakly depends on A^^ and kp- Applying the same ratio to ^ at the 
thermodynamic limit, we then obtain 

^(r/ ^ -0.8) = 0.728(37), 

i{r] ^ -0.4) = 0.546(34), (13) 
^{ri ^ -0.2) = 0.428(30). 

The three values of ^ of Eq. (fT3|l are best fitted by 

^ = 0.292(24) - 0.795(33) r] - 0.271(21) r^^, (14) 

which yields ^ = 0.292(24) at the unitary limit by setting 77 — 0. Note that Eq. f|T^ is 
similar to ^ = 0.306(1) -0.805(2)r/- 0.63(3)7^2 ^ and gives d^/d7]\r,=o = -0.795(33) 

close to -1.0(1) in Ref. Q. 

Figure [3] shows the three ^ values of Eq. (fT3|) and their best fit Eq. (IT^ , together with 
^ at the unitary limit. In the figure, we also show several results of the 1] dependence of ^ 
reported in the literature. They include the two types of calculations: (i) analytical and (ii) 
numerical. 



(i) This type is an e-expansion calculation at the next-to-leading order (NLO) 231] about 



four dimension (shown by the dash-dotted curve) . Note that a recent next-to- next-to-leading 



order (NNLO) calculation at 77 = 2J] shows the appearance of a problematic In e contri- 
bution, but it is argued to be infrared manageable. 

(ii) The Monte Carlo types of calculation include the lattice Monte Carlo calculation with 



the symmetric heavy-light ansatz 22] (shown by the dotted curve), the diffusion Monte Carlo 



method 251] (cross symbols), the fixed- node Green's function Monte Carlo method at 77 < 



(solid up-triangles) [26], and at 77 = (sohd down-triangle) [27|], and the determinantal 



quantum Monte Carlo method (solid squares) 28|, |29|]. Note that the previously reported ^ 



by the fixed-node Green's function Monte Carlo calculations [30|, l3ll (solid diamonds) were 
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somewhat larger, being near that of Ref . 25|] . We further note that no thermodynamic hmit 
is taken in these Monte Carlo calculations. 




FIG. 3: 1] = l/^kpao) dependence of the energy parameter ^. The solid circles and the dashed curve 
are our Monte Carlo data and the best-fit, Eq. (jl4p . respectively. Our ^ extrapolated to the unitary 
limit is shown by the open circle. For comparison, ^ obtaine d by other works are also shown: the 
dash-dotted curve is by the next-to-leading order e expansion 
Monte Carlo calculation with the symmetric heavy-light ansatz 
fixed-node diffusion Monte Carlo calculation 



23] , the dotted curve is by the lattice 



22| . the cross symbols are by the 



25l | , the solid up- and down-triang. 



node Green's function Monte Carlo calculations at < [26] and at r/ = [271]. respectively, and 



the solid squares are by the determinantal quantum Monte Carlo calculation 



es are by the fixed- 



28|, 



29|. 



In Fig.m we compare our result with various values of the unitary-limit ^ in the literature. 
The figure is made by expanding a similar figure (Fig. 14) in Ref. |24]. On the left-hand 
side of the figure, the values of ^ determined by atomic Fermi-gas experiments are shown 
by solid circles: ^ = 0.74(7) fl. ^ = 0.34(15) 33|, ^ = 0.32l°:}i^ [34], ^ = 0.46(5) 35], and 
e = 0.51(4) [36]. 

In the middle of Fig. HI the values of C, obtained by various Monte Carlo calculations 
are shown^by the solid sguares: ^ < 0.42(1) [27], ^ = 0_^(3) 0, 0.07 < ^ ^0^2 11|, 



^ ^ 0.44 



37|, ^ ^ 0.28 [22], ^ = 0.292(12) and 0.329(5) [12], and ^ = 0.37(5) 



281, 



29(. For 



comparison, our result, Eq. (fT4|) . is shown by the diamond. 
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FIG. 4: ^ at the unitary limit. The values of ^ reported in the literature are shown in three 
groups from left to right: those determined by atomic Fermi-gas experiments (solid circles); by 
various Monte Carlo calculations (solid squares); and by e expansions (triangles). The scale of the 
horizontal axis has no significance but for separating data, .^'s by the e-expansion method on the 
right-hand side are divided into three subgroups as explained in the text. Our ^ from Eq. (jl4p 
is shown by the diamond at the right-most location in the Monte Carlo group. See the text for 
the reference to each value of ^. The figure is an expanded version of a similar figure (Fig. 14) in 



Ref. 



On the right-hand side of the figure, the values of ^ by the e expansions are shown by 
triangles in three groups: from left to right, (1) the NLO e expansion [38], (2) the Borel-Pade 
approximation between the NLO expansions about four and two dimensions 39||, (3) the 
Borel-Pade approximation between the NNLO expansion around four dimensions and the 
NLO expansion about two dimensions |2^, and (4) the Borel-Pade approximation between 



the NNLO expansions about four and two dimensions 



4Q| 



Figure m shows that our value of ^, 0.292(24), is relatively small among the values shown. 



IV. PAIRING GAP A 

The pairing gap A at the unitary limit may be simply related to Eg,s. as ~ 2Eg,_ 
To examine the relationship, we determined A for rj ^ hj extrapolation. 



3Q| 
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In Fig. [5l we show A in the thermodynamic and continuum hmits for the three values of 

n 

7] by taking from our previous work at the LO p^]. The figure also shows the best-fit curve 



to the three rj values using the quadratic function 

A 



0.384(30) + 0.303(27) r/ + 0.046(37) r]\ 



(15) 



from which A/e^? = 0.384(30) is determined by setting 77 = 0. The figure includes other 
quantum Monte Carlo results: the solid squares are by Bulgac et al. [28|, l29|], the solid up- 
triangles are by Gezerlis et al. [26], the solid down-triangle is by Carlson et al. 27 1, and the 



solid diamonds are by Chang et al. 



to Fig. 14 of Ref. 
Refs. 28 



3l| . Figure [5] is drawn similar to Fig. 1 of Ref. 



28 



and 



291]. Note that the data by Chang et al. at 77 < are those quoted in 



The relation between A and Eg s. is found from that of A and r] by eliminating ep from 
Eqs. (11 01) and (|T5l) and by using of Eq. at r/ = 0. Table II summarizes our result. Our 
A/Eg,s. at the unitary limit is 2.19(35) and roughly confirms ~ 2 as suggested in Ref. 30] . 
Note that other quantum Monte Carlo calculations in the literature also yield similar values: 
A/Eg,s. = 2.5(3) 



28 



29|, 2.3(1) |3l|, and 2.0(2) |27|. 



TABLE II: A/Eg,s. for various values of r/ 



A/E, 



g.s. 



-0.7997 0.172(9) 0.728(37) 0.391(40) 

-0.3999 0.271(16) 0.546(34) 0.827(99) 

-0.1999 0.326(22) 0.428(30) 1.27(18) 

0.0000 0.384(30) 0.292(24) 2.19(35) 



V. CRITICAL TEMPERATURE 



In the previous two sections, we examined the physical quantities at zero temperature. 
Thermodynamics at the unitary limit is of much interest [7, 



28 



29l |. In this section, we 



examine a representative thermodynamic quantity: the critical temperature of the phase 
transition, T^, at the unitary limit. 
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< 




FIG. 5: Pairing gap A in the unit of ep as a function of rj = l/{kFao). The sohd circles are our 
Monte Carlo data, shown with statistical uncertainties The dashed curve is the best fit to the 
data by use of a quadratic function of r/. A at the unitary limit r] ^ determined by extrapolation 
is shown by an open circle. The solid squares, the solid up-triangles, the solid down-triangle, and 
the solid diamonds are the quantum Monte Carlo data by Bulgac et al. 



28 



al. 



2a|, by Carlson et al. [27], and by Ch ang et al. 



et al. at ?7 < are those quoted in Refs. 



29| 



29l |. by Gezerlis et 



3ll |. respectively. Note that the data by Chang 



Figure [6] shows Tc in the thermodynamic and continuum limits for the three values of r] 
by taking from our previous work at the LO The best fit to the data is also shown in 
the figure by the dashed curve. It is given by 

T 

— = 0.189(12) + 0.149(22) r] + 0.069(36) t]^, (16) 
ep 

which yields the extrapolated value Tc/ep = 0.189(12) aX rj = 0, or at the unitary limit. 
Jot comparison, in the figure we show Tc/ep by a recent quantum Monte Carlo calculation 
29 1 . While our values are somewhat larger than those of Ref. j^], both show similar rj 
dependence on T^. This rj-Tc dependence is also observed by an e-expansion calculation j41|. 
Our Tc is tabulated for various values of t] in Table III. 

In Fig. [71 we summarize various values of at the unitary limit that are reported in the 
literature, together with ours for comparison. The figure includes the two groups of the Tc 
determination: by Monte Carlo calculations and by the e-expansion method. 
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FIG. 6: r] dependence of the critical temperature Tc, shown in the unit of ep. Our Monte Carlo 
data at the thermodynamic and continuum limits are shown by solid circles, and the extrapolated 
unitary-limit point is shown by an open circle. The extrapolation is made using the fit function, 
Eq. ()16p . For comparison, the solid squares show a recent quantum Monte Carlo result 29 1. 



TABLE III: Tc/ep for various values of r]. 

V Tc/ep 
-0.7997 0.114(23) 
-0.3999 0.141(10) 
-0.1999 0.162(11) 
0.0000 0.189(12) 



The Monte Carlo results include T^/ep = 0.23(2) by the determinantal quantum Monte 



Carlo method on A^, = 6'^ and 8^ lattices 



37 



42l |: Tc/ep < 0.15(1) by the determinantal 



quantum Monte Carlo method on A^^ = 6^, 8^, and 10^ lattices 29j; T^/ep < 0.14 by the 
hybrid Monte Carlo method on A^, = 4^, 5^, and 6^ lattices [111]; Tjep = 0.152(7) by a 
diagrammatic determinantal quantum Monte Carlo method with the finite-size scaling {g]; 
and Tc/ep ~ 0.25 by a restricted path integral Monte Carlo method 43|. 

The e-expansion results \^ include T^/ep ~ 0.249 and T^/ep ~ 0.153 by up to the NLO 
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FIG. 7: Critical temperature Tc at the unitary limit appearing in the literature and our extrapolated 



Tc (all shown in the unit of ep)- The solid squares (with statistical 



side of the figure are by Monte Carlo calculations a, 111 , 



291, 



37 



uncertainties) on the left-hand 



43], and the solid triangles 



411 ] ■ Our Tc is shown by the solid diamond 



on the right-hand side are by e expansion calculations 
with the statistical uncertainty. The horizontal scale is used merely for the separation of each Tc 
value. The Tc's by the e expansion calculations are grouped into two in the same way as in Fig. [H 
See the text for the reference of each Tc. 



in the e expansion about the four and two dimensions, respectively; and Tc/ep = 0.173 and 
0.192 by the Borel-Pade approximation between the four and two dimensions. As noted in 
the case of the ^ determination, the e expansion seems to indicate a possibly problematic 
behavior at higher orders 24]. 

Our Tc at the thermodynamic and continuum limit is also shown in the Monte Carlo group 
of Fig. [71 Our result Tc/ep ~ 0.189(12) seems to be consistent with others, especially with 
that of Ref. 6|], which is the only other Tc obtained in the thermodynamic and continuum 
limits. 



VI. CONCLUSION 

We extrapolate our Monte Carlo results for low-density neutron matter jl| to the unitary 
limit for the following quantities: the ground-state energy Eg s., the pairing gap A at T ^ 0, 
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and the critical temperature of the normal-to-superfluid phase transition Tc. All quantities 
show a smooth extrapolation to the limit in terms of 77 = l/{kFao)- Although no accurate 
determination of these quantities is yet available, our extrapolated values are in reasonable 
agreement with those in the literature. Our successful extrapolation suggests that much of 
the physics of low-density neutron matter [of about (10^'^ — 10^^)po] is similar to the physics 
of the attractively interacting fermion system at the unitary limit. 
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